# This file contains the code to reproduce graphs, tables, and quantitative statements from section 3 of the article ``Environmental Policy and Directed Technological Change: Evidence from the European carbon market''. For convenience, the script is organised as far as possible to follow the order in which results are presented in the article.

############### Preamble ###############
require(strucchange) # If you do not have this package installed, run the command "install.packages("strucchange", dependencies=T)" first.
rm(list = ls())
load("EUETS_innovation_public.Rdata")

############### Replication code for section 3 ###############
## Code to reproduce graphs, tables, and quantitative claims in section 3:
# "The EPO was set up in 1978. Since then, over 2.5 million patents have been filed with the EPO, of which just over 50,000 (or 2%) have been classified as low-carbon inventions."
sum(descriptive$tot_epo_pat) # Total number of patents filed at the EPO
sum(descriptive$tot_epo_green_pat_Y02) # Total number of low-carbon patents filed at the EPO
sum(descriptive$tot_epo_green_pat_Y02)/sum(descriptive$tot_epo_pat) # Share of low-carbon patents filed at the EPO

# "the firms in our data set account for just over 95% of all patents filed at the EPO"
sum(descriptive$total_epo_patents_orbis4)/sum(descriptive$tot_epo_pat) # Share of EPO patents successfully linked to BvD ID numbers in our database

# Figure 1
quartz(width=7, height=4.5, "share7809", file = "share7809.pdf", type = "pdf")
plot(descriptive$year, (descriptive$tot_epo_green_pat_Y02 / descriptive$tot_epo_pat)*100, type="l", xlab="Year", ylim=c(0,5), ylab="Share of all patents (in %)")
abline(v=c(2005),lty=3)
text(x=2005,y=4.9, pos=2, cex=1, "EU ETS")
text(x=1981,y=2.6, "Low-carbon")
lines(descriptive$year, (descriptive$tot_epo_pct_pat / descriptive$tot_epo_pat)*100, lty = 5)
text(x=1988,y=0.3, "Pollution control")
dev.off()

# "A simple Chow test strongly rejects the hypothesis that there is no structural break in 2005 (P < 0.001)."
share <- ts(data=(descriptive$tot_epo_green_pat_Y02 / descriptive$tot_epo_pat), start=1978)
fs <- Fstats(share ~ 1)
breakdates(breakpoints(fs)) # 2005 is the top candidate for a structural break
fs <- Fstats(share ~ 1, from=28, to=28)
sctest(fs) # P-value associated with 2005

# "this pattern is robust to using an expanded definition of "low-carbon technologies""
share <- ts(data=(descriptive$tot_epo_green_pat_Y02plus / descriptive$tot_epo_pat), start=1978)
fs <- Fstats(share ~ 1)
breakdates(breakpoints(fs)) # 2005 is still the top candidate for a structural break
fs <- Fstats(share ~ 1, from=28, to=28)
sctest(fs) # P-value associated with 2005

# "the share of patents protecting non-greenhouse gas “pollution control technologies”, as defined by Popp (2006), which does not display the same structural break (one cannot reject the hypothesis of no structural break in 2005 at conventional significance levels)."
share <- ts(data=(descriptive$tot_epo_pct_pat / descriptive$tot_epo_pat), start=1978)
fs <- Fstats(share ~ 1)
breakdates(breakpoints(fs)) # Top candidate for structural break no longer 2005
fs <- Fstats(share ~ 1, from=28, to=28)
sctest(fs) # P-value associated with 2005

# Figure 2
quartz(width=7, height=4.5, "oil", file = "oil.pdf", type = "pdf")
par(mar=c(5,4,4,5)+.1)
plot(descriptive$year, (descriptive$tot_epo_green_pat_Y02 / descriptive$tot_epo_pat)*100,type="l", ylim=c(0,5), xlab="", ylab="Share of low-carbon patents (in %)")
text(x=1986,y=0.8, "Low-carbon patents")
text(x=1986,y=4, "Crude oil price")
text(x=2005,y=4.9, pos=2, cex=1, "EU ETS")
abline(v=c(2005),lty=3)
par(new=TRUE)
plot(descriptive$year, descriptive$oil_price_2010_usd,type="l",xaxt="n", ylim=c(0,100), yaxt="n",lty = 5,xlab="Year",ylab="", frame.plot =FALSE)
axis(4)
mtext("2010 USD",side=4,line=3)
dev.off()

# Table 1
table1 <- data.frame(country=c("AT", "BE", "CZ", "DK", "EE", "FI", "FR", "DE", "IE", "LT", "LU", "NL", "PL", "PT", "SK", "ES", "SE", "UK"), stringsAsFactors=F)
for (i in 1:nrow(table1)) {
	table1$numberofinstallations[i] <- nrow(ETSinstallations[ETSinstallations$registry_code==table1$country[i],]) # Count the number of EU ETS installations
}
for (i in 1:nrow(table1)) {
	table1$mtonnesofemissions[i] <- (sum(ETSinstallations$allocated_2005[ETSinstallations$registry_code==table1$country[i]], na.rm=T) + sum(ETSinstallations$allocated_2006[ETSinstallations$registry_code==table1$country[i]], na.rm=T) + sum(ETSinstallations$allocated_2007[ETSinstallations$registry_code==table1$country[i]], na.rm=T))/(10^6) # Calculate Mtonnes of allocated emissions
}
for (i in 1:nrow(table1)) {
	table1$percentofinstallationscovered[i] <- 100*nrow(ETSinstallations[ETSinstallations$registry_code==table1$country[i] & ETSinstallations$matched==1,]) / nrow(ETSinstallations[ETSinstallations$registry_code==table1$country[i],]) # Calculate percent of installations operated by identified firms
}
for (i in 1:nrow(table1)) {
	table1$percentofemissionscovered[i] <- 100*(sum(ETSinstallations$allocated_2005[ETSinstallations$registry_code==table1$country[i] & ETSinstallations$matched==1], na.rm=T) + sum(ETSinstallations$allocated_2006[ETSinstallations$registry_code==table1$country[i] & ETSinstallations$matched==1], na.rm=T) + sum(ETSinstallations$allocated_2007[ETSinstallations$registry_code==table1$country[i] & ETSinstallations$matched==1], na.rm=T)) / (sum(ETSinstallations$allocated_2005[ETSinstallations$registry_code==table1$country[i]], na.rm=T) + sum(ETSinstallations$allocated_2006[ETSinstallations$registry_code==table1$country[i]], na.rm=T) + sum(ETSinstallations$allocated_2007[ETSinstallations$registry_code==table1$country[i]], na.rm=T)) # Calculate percentage of regulated emissions from installations operated by identified firms
}
# Add totals
table1 <- rbind(table1, c("Total", nrow(ETSinstallations[ETSinstallations$registry_code %in% table1$country,]), (as.numeric(sum(ETSinstallations$allocated_2005[ETSinstallations$registry_code %in% table1$country], na.rm=T)) +  as.numeric(sum(ETSinstallations$allocated_2006[ETSinstallations$registry_code %in% table1$country], na.rm=T)) + as.numeric(sum(ETSinstallations$allocated_2007[ETSinstallations$registry_code %in% table1$country], na.rm=T))) / (10^6), 100*nrow(ETSinstallations[ETSinstallations$registry_code %in% table1$country & ETSinstallations$matched==1,]) / nrow(ETSinstallations[ETSinstallations$registry_code %in% table1$country,]), 100*(as.numeric(sum(ETSinstallations$allocated_2005[ETSinstallations$registry_code %in% table1$country & ETSinstallations$matched==1], na.rm=T)) +  as.numeric(sum(ETSinstallations$allocated_2006[ETSinstallations$registry_code %in% table1$country & ETSinstallations$matched==1], na.rm=T)) + as.numeric(sum(ETSinstallations$allocated_2007[ETSinstallations$registry_code %in% table1$country & ETSinstallations$matched==1], na.rm=T))) / (as.numeric(sum(ETSinstallations$allocated_2005[ETSinstallations$registry_code %in% table1$country], na.rm=T)) +  as.numeric(sum(ETSinstallations$allocated_2006[ETSinstallations$registry_code %in% table1$country], na.rm=T)) + as.numeric(sum(ETSinstallations$allocated_2007[ETSinstallations$registry_code %in% table1$country], na.rm=T)))))
table1 <- rbind(table1, c("Total EU ETS", nrow(ETSinstallations), (as.numeric(sum(ETSinstallations$allocated_2005, na.rm=T)) +  as.numeric(sum(ETSinstallations$allocated_2006, na.rm=T)) + sum(as.numeric(ETSinstallations$allocated_2007, na.rm=T))) / (10^6), 100*nrow(ETSinstallations[ETSinstallations$registry_code %in% table1$country & ETSinstallations$matched==1,]) / nrow(ETSinstallations), 100*(as.numeric(sum(ETSinstallations$allocated_2005[ETSinstallations$registry_code %in% table1$country & ETSinstallations$matched==1], na.rm=T)) +  as.numeric(sum(ETSinstallations$allocated_2006[ETSinstallations$registry_code %in% table1$country & ETSinstallations$matched==1], na.rm=T)) + as.numeric(sum(ETSinstallations$allocated_2007[ETSinstallations$registry_code %in% table1$country & ETSinstallations$matched==1], na.rm=T))) / (as.numeric(sum(ETSinstallations$allocated_2005, na.rm=T)) + as.numeric(sum(ETSinstallations$allocated_2006, na.rm=T)) + sum(as.numeric(ETSinstallations$allocated_2007, na.rm=T)))))

# Figure 3
quartz(width=7, height=4.5, "relativeshares", file = "relativeshares.pdf", type = "pdf")
plot(descriptive$year, (descriptive$greenpat_epo_ets / descriptive$total_patents_ets) *100, type="l", xlab="Year", ylim=c(0,5), ylab="Share of low-carbon patents (in %)")
text(x=1985,y=3.5, "EU ETS firms")
lines(descriptive$year, (descriptive$greenpat_epo_nonets / descriptive$total_patents_nonets)*100, lty = 5)
abline(v=c(2005),lty=3)
text(x=1984,y=0.6, "non-EU ETS firms")
text(x=2005,y=4.9, pos=2, cex=1, "EU ETS")
dev.off()

# Footnote: "the average number of citations received by low-carbon patents filed by EU ETS companies since 2005 does not significantly differ from those filed by non-EU ETS companies."
plot(descriptive$year[descriptive$year>2004], (descriptive$greenpat_epo_cit_ets[descriptive$year>2004] / descriptive$greenpat_epo_ets[descriptive$year>2004]), type="l", xlab="Year", ylab="Citations per low-carbon patent", ylim=c(0,1))
# text(x=1985,y=3.5, "EU ETS firms")
lines(descriptive$year[descriptive$year>2004], (descriptive$greenpat_epo_cit_nonets[descriptive$year>2004] / descriptive$greenpat_epo_nonets[descriptive$year>2004]), lty = 5)
text(x=2007, y=0.2, "non-EU ETS firms")
text(x=2008, y=0.4, pos=2, cex=1, "EU ETS")

# Footnote: "Similarly, the size of low-carbon patent families is the same for EU ETS and non-EU ETS companies."
plot(descriptive$year[descriptive$year>2004], (descriptive$greenpat_epo_fam_ets[descriptive$year>2004] / descriptive$greenpat_epo_ets[descriptive$year>2004]), type="l", xlab="Year", ylab="Average family size of low-carbon patent", ylim=c(0,6))
# text(x=1985,y=3.5, "EU ETS firms")
lines(descriptive$year[descriptive$year>2004], (descriptive$greenpat_epo_fam_nonets[descriptive$year>2004] / descriptive$greenpat_epo_nonets[descriptive$year>2004]), lty = 5)
text(x=2008, y=5, "non-EU ETS firms")
text(x=2007, y=4, pos=2, cex=1, "EU ETS")

# "EU ETS firms filed 2,189 low-carbon patents in 2005–2009, compared to 972 patents in the 5 preceding years (an increase of 125%)..."
sum(descriptive$greenpat_epo_ets[descriptive$year>2004]) # Low-carbon patents 2005–2009
sum(descriptive$greenpat_epo_ets[descriptive$year>1999 & descriptive$year<2005]) # Low-carbon patents 2000-2004
(((sum(descriptive$greenpat_epo_ets[descriptive$year>2004]) - sum(descriptive$greenpat_epo_ets[descriptive$year>1999 & descriptive$year<2005])) / sum(descriptive$greenpat_epo_ets[descriptive$year>1999 & descriptive$year<2005])))*100 # Percentage change

# "... while non-EU ETS firms filed 19,841 and 12,037 low-carbon patents in the corresponding periods (an increase of 65%)."
sum(descriptive$greenpat_epo_nonets[descriptive$year>2004]) # Low-carbon patents 2005–2009
sum(descriptive$greenpat_epo_nonets[descriptive$year>1999 & descriptive$year<2005]) # Low-carbon patents 2000-2004
(((sum(descriptive$greenpat_epo_nonets[descriptive$year>2004]) - sum(descriptive$greenpat_epo_nonets[descriptive$year>1999 & descriptive$year<2005])) / sum(descriptive$greenpat_epo_nonets[descriptive$year>1999 & descriptive$year<2005])))*100 # Percentage change

# "we can naively estimate how many low-carbon patents the EU ETS has added so far: 2,189 - 1.65 × 972 = 585.2."
sum(descriptive$greenpat_epo_ets[descriptive$year>2004]) - (sum(descriptive$greenpat_epo_nonets[descriptive$year>2004]) / sum(descriptive$greenpat_epo_nonets[descriptive$year>1999 & descriptive$year<2005]))*sum(descriptive$greenpat_epo_ets[descriptive$year>1999 & descriptive$year<2005]) # Number in text is approximate, and includes minor rounding error as a result of assuming a 65% increase instead of a 64.8334% increase.

# "This amounts to a 2.6% increase in the number of low-carbon patents at the EPO compared to what it would have been without the EU ETS."
100*(sum(descriptive$greenpat_epo_ets[descriptive$year>2004]) - (sum(descriptive$greenpat_epo_nonets[descriptive$year>2004]) / sum(descriptive$greenpat_epo_nonets[descriptive$year>1999 & descriptive$year<2005]))*sum(descriptive$greenpat_epo_ets[descriptive$year>1999 & descriptive$year<2005])) / (sum(descriptive$tot_epo_green_pat_Y02[descriptive$year>2004]) - (sum(descriptive$greenpat_epo_ets[descriptive$year>2004]) - (sum(descriptive$greenpat_epo_nonets[descriptive$year>2004]) / sum(descriptive$greenpat_epo_nonets[descriptive$year>1999 & descriptive$year<2005]))*sum(descriptive$greenpat_epo_ets[descriptive$year>1999 & descriptive$year<2005]))) # Number in text is approximate, and includes minor rounding error as a result of assuming a 65% increase instead of a 64.8334% increase.

# "while only 1 in about 5,500 firms is EU ETS regulated, they account for roughly 1 in 12 low-carbon patents filed in the 5 years before the EU ETS launched."
30000000 / 5568 # Ratio of EU ETS to non-EU ETS firms in text is apprioximate.
sum(descriptive$tot_epo_green_pat_Y02[descriptive$year>1999 & descriptive$year<2005]) / sum(descriptive$greenpat_epo_ets[descriptive$year>1999 & descriptive$year<2005]) # Ratio of EU ETS to non-EU ETS low-carbon patents in text is approximate, and includes minor rounding error.